Revealing Interface Polarization Effects on the Electrical Double Layer with Efficient Open Boundary Simulations under Potential Control

A major challenge in modeling interfacial processes in electrochemical (EC) devices is performing simulations at constant potential. This requires an open-boundary description of the electrons, so that they can enter and leave the computational cell. To enable realistic modeling of EC processes under potential control we have interfaced density functional theory with the hairy probe method in the weak coupling limit (Zauchner et al. Phys. Rev. B2018, 97, 045116). Our implementation was systematically tested using simple parallel-plate capacitor models with pristine surfaces and a single layer of adsorbed water molecules. Remarkably, our code’s efficiency is comparable with a standard DFT calculation. We reveal that local field effects at the electrical double layer induced by the change of applied potential can significantly affect the energies of chemical steps in heterogeneous electrocatalysis. Our results demonstrate the importance of an explicit modeling of the applied potential in a simulation and provide an efficient tool to control this critical parameter.

3) To better elucidate water polarization at charged surfaces, considering a water monomer model might offer a more straightforward analysis.4) In presenting calculated capacitance values in Table 1, adopting the unit "uF/cm2" would facilitate easier comparison with both experimental and theoretical studies in the field.
5) The discrepancy in slab separation distances-10 angstroms for some models and 20 angstroms for others-within the method section warrants explanation.Clarifying the rationale behind these choices would aid in understanding the modeling approach.

Comments to the Author
The manuscript by Buraschi, Horsfield and Cucinotta discusses the implementation of the hairy probes formalism, in which each atom at an interface of interest is connected to a lead imposing a desired potential, into the CP2K package.This approach is then applied to study properties of a water bi-layer/Pt(111) interface.The work is interesting, but there are limitations I elaborate on below, which make the manuscript unsuitable for publication in its present form.After these have been addressed, the manuscript should be re-evaluated.

1)
The authors write that the leads are typically coupled only to atoms in the contact region.Obviously, this applies to the electrode and is needed to achieve a potential drop at the interface.They also mention the implementation of "solutions probes" needed to avoid a "nonphysical situation where molecules in solution have no electrons".Do I understand correctly that this "solution leads" couple to water molecules, so that this coupling leads to a surplus/deficiency of electrons on the water molecules?What are the charges on the water molecules?

2)
Some of the applied potentials are rather high and I would expect reactions to occur.Yet the shown effects are molecule desorption or reorientation.It is not clear to me, whether reactions, for example water dissociation, can occur.How will these be affected by the "solution leads"?

3)
If I understand the set-up correctly, the investigated interface consists of a Pt(111) slab in contact with two water layers, which in turn are connected to a vacuum region.What about a bulk water region?The properties of an adsorbed water bi-layer in vacuum are not necessarily the same as the properties of an adsorbed water bi-layer in contact with water.Yet the shown results are often compared to calculations in which an extended water region is present.Is such a comparison applicable?Are the comparisons of timing of calculations using an extended bulk water region fare?Obviously any calculations which does not have to take care of equilibrating a bulk water region will be faster.

4)
I find the discussion of "The response of a water bilayer to the charge of electrode potential" to be a bit selective on the comparison to available literature and selection of which articles to cite.The understanding of electrochemical interfaces has evolved since the cited 2004 article and people are aware that changes in the electrode potential affect the polarisation, orientation and adsorption/desorption behaviour of water molecules at the interface, even if electron transfer reactions have not yet occurred.

5)
In the discussion of the set-up the authors mentioned two numbers for the separation between the Pt slabs -10 Ang and 20 Ang.It is not clear to me, whether the separation of 20 Ang applies to all the models containing water or only to the ones with a reduced number of water molecules.It is also not clear why it was necessary to change the separation length.

6)
Why is an entropy term included only for the gas-phase molecule?What about the water molecules at the interface?What about entropic contribution by the bi-layer water molecules?Surely this will be lower, but not necessarily negligible.Do the authors have an estimate about their magnitude and whether neglecting them is justified?

7)
I have difficulties with the formulation "the potential increases" (p.11).I assume the authors mean "becomes more positive/negative"?

8)
What does a "full coverage water bilayer" mean -a 12 water molecules structure in the first and 12 water molecules in the second layer?What is the meaning of "low water coverage structures"?Please be more precise.

9)
I believe that a citation is missing on p.12 when the authors say that their "results also align with recent literature".

10)
What is the meaning of E_h?

11)
Please be more specific when writing that the "HP-DFT formalism in calculations did not exhibit a significant deceleration of the SCF cycle."Some numbers were provided in the SI comparing timings for one cycle.What is not clear to me is, how long does the overall equilibration of the system take.After all, one needs to equilibrate the various leads.How efficient the overall performance is, would likely very much depend on the implementation and used algorithm.

12)
While I can guess which DOS in Fig. 2b relates to the left and which to the right lead, it would be helpful to also write this explicitly in the figure.

13)
It would be helpful to incorporate the labels of the graphs in Fig. 4b into the figure.

14)
How are distances measured?Do the authors use the centre of mass of a water molecule or maybe the oxygen in a water molecule?
Author's Response to Peer Review Comments: Dear Editor, Thank you for considering our manuscript for publica:on in "The Journal of Physical Chemistry Le@ers" and for forwarding the comments of the referees.We have carefully reviewed their comments and sugges:ons, and we are pleased to submit our detailed responses addressing all technical and analy:cal issues raised.
A@ached, you will find our point-by-point response to the referees' comments (highlighted in yellow), along with a list of the addi:onal changes made to the manuscript in light of their feedback, with addi:ons marked in blue for the main text and in red for the supplementary material.
We have included two versions of the revised Manuscript and Supplementary Informa:on: The documents " main_reviewed.pdf" and " Suppor:ng_reviewed.pdf" highlight the changes made compared to the previous version of the manuscript.They are uploaded as part of the manuscript zip file.Please note that references to line numbers, tables, figures, and equa:ons correspond to these revised documents.
Addi:onally, we have provided polished versions of both the manuscript and the Supplementary Informa:on, separately.
A list of addi:onal formaUng modifica:ons is reported at the end of this document.
We believe that our revisions address the concerns raised by the referees and enhance the overall quality of our manuscript.We hope you will now find our work suitable for publica:on in "The Journal of Physical Chemistry Le@ers."Yours, Sincerely, Margherita Buraschi, Andrew Horsfield and Clo:lde Cucino@a Reviewer #1 1 The manuscript did not clearly outline the advantages of the HP DFT algorithm compared with other constant poten:al algorithms.
The HP-DFT algorithm exhibits its key advantage in its efficiency.The Hairy-Probes formalism was easily implemented in the DFT code by modifying only the way the occupa:on numbers are computed.Given that this modifica:on cons:tutes the sole devia:on from the regular SCF cycle, most of the required variables are already computed and provided by the rest of the DFT code.Consequently, the :ming of an SCF cycle in an HP-DFT calcula:on does not showcase any significant slowdown with respect to the :ming of an SCF cycle in standard DFT calcula:on.This is highlighted in the SI, in the code efficiency sec:on.
To give a more quan:ta:ve analysis of the efficiency, the SI has been modified as follows: • lines 8 to 16 in Suppor0ng_reviewed.pdf: "For example, a model such as the Pt(111)(6×6 Another feature which makes HP-DFT a one-of-a-kind formalism is that it allows for a mul:-terminal setup.Although this is not showcased in the manuscript, poten:ally any number of probe sets with different electrochemical poten:als can be used.This paper served mostly to benchmark the formalism, so we worked on simple systems.Furthermore, the way the HP-DFT has been implemented, the cell is charge neutral and the number of par:cles does not vary, ensuring that the energy is well defined.Finally, the probes directly control the electron electrochemical poten:als, thereby providing a proper descrip:on of the electrons in the system.
To make this point clearer, we have modified the manuscript as follows: • lines 67 and 68 of main_reviewed.pdf:"HP-DFT proved to be lightweight and highly efficient, showcasing computaGonal costs comparable to that of standard DFT calculaGons.".• lines 162 to 170 of main_reviewed.pdf:"Because of the way the HP-DFT was implemented, the cell is charge neutral.This ensures that the energy is well defined.Finally, the probes impose well-defined electron electrochemical potenGals, thereby providing a proper descripGon of the electrons in the system.Finally, another feature which makes HP-DFT a one-of-a-kind formalism is its capability to accommodate a mulG-terminal setup.While exploring this aspect is beyond the scope of this manuscript, any number of probe sets with different electrochemical potenGals could potenGally be employed.As the focus of this paper primarily centred on benchmarking the formalism, our efforts were directed towards studying systems with a two-terminal setup for simplicity".• lines 173 to 180 of main_reviewed.pdf: and "Overall, the key feature of HP-DFT is its efficiency.The integraGon of the Hairy-Probes formalism into the DFT code was easily achieved, requiring only a modificaGon in the way the occupaGon numbers are calculated.As this is the sole deviaGon from the standard SCF cycle, the majority of the necessary variables are readily computed and supplied by the exisGng DFT code infrastructure.Consequently, the use of the HP-DFT formalism in calculaGons did not exhibit a significant deceleraGon of the SCF cycle in comparison to convenGonal DFT algorithms.Further informaGon regarding the code's efficiency is available in Figure S1 of the SupporGng Material." .

(a)
In the "HP-DFT formalism" sec:on, the authors men:on a coupling strength parameter Γ !, but this parameter is not observed in subsequent formulas.(b) It has also not been clarified how to calculate the average electron electrochemical poten:al ̅ .(c) Addi:onally, the authors applied  " and  # based on the assump:on of ∆ " = −∆ # .However, this assump:on may not be appropriate, as the simula:on model is not symmetrical and water molecules are present only on one side of the electrode.
(a) The HP-DFT formalism is implemented in the limit of weakly coupled probes.In the original formula:on of strongly coupled probes, Γ !appears both in the formula for the occupa:on numbers and in the Hamiltonian of the system.In the weak coupling limit, it is assumed that Γ !→ 0, so this parameter does not appear in our formulas.
To clarify this point, the manuscript has been modified as follows: lines 81 to 85 of main_reviewed.pdf:"In the limiGng case of probes weakly coupled to the system, Γ p → 0, thus it does not appear in the central formulas.It was shown that this weak coupling limit is suitable for describing EC systems as the charge is carried from one electrode to the other by ions in the electrolyte, whose diffusion rate sets the electron conducGon rate between electrodes.".
(b) The average Fermi level, ̅ , is calculated as the average among the local Fermi levels induced by the probes' EC poten:al.̅ is the level which ensures charge neutrality in the system.Essen:ally, we set the value of the EC poten:al difference, ∆ = |∆ " − ∆ # |, and define ̅ as the midpoint between the two Fermi levels.This value is then adjusted to enforce charge neutrality effec:vely.
To clarify this point, the manuscript has been modified as follows: lines 113 to 115 of main_reviewed.pdf:"The level ̅ is calculated as the average among the local Fermi levels induced by the probes' EC potenGal; as such, this is the level which ensures charge neutrality in the system.".
(c) It is correct that the simula:on model is not symmetrical and, therefore, the electrodes poten:als would not, on their own, sa:sfy this rela:on.However, ∆ " and ∆ # are defined as the probes' electrochemical poten:als referred to ̅ and they are set by the user.In this case, therefore, it is not assumed that ∆ " = −∆ # ; these are the values that have been chosen to assign them.In principle we have three electrochemical poten:als values to set.We currently only use two criteria to set them (charge neutrality and applied bias), so we need an addi:onal rule (as described above) to uniquely define all three values.We note that  is also used to set the EC poten:al for the solu:on probes.
To clarify this point, the manuscript has been modified as follows: line 159 of main_reviewed.pdf:"were chosen to saGsfy" 3 The interface model used in this manuscript contains a small number of water molecules only at one side electrode, while different electrochemical poten:als are applied to both electrodes.In this model, water molecules may be unable to screen the excess charges on both electrodes, which could lead to serious issues with periodic boundary condi:ons and result in unreliable simula:on results.To address this issue, the author should consider an implicit solvent model or dipole correc:on, or widen the electrode spacing, incorpora:ng a sufficient number of water molecules, and conduc:ng AIMD calcula:ons.
While it is en:rely true that there is not enough water to screen the dipole introduced by the charged slabs, this does not represent an issue for our poisson solver as periodic boundaries are not used in the direc:on perpendicular to the slabs' surfaces and parallel to the dipole.Instead, an implicit Posson solver is employed.We acknowledge that this may not come across clearly in the text and we have modified the manuscript as follows: lines 151 to 154 of main_reviewed.pdf:"The introducGon of an EC potenGal difference induces a dipole within the system.For this reason, periodic boundary condiGons were applied only in the direcGons parallel to the surfaces of the plates (x and y), while in the direcGon perpendicular to them (z) an implicit Poisson solver was used.
Due to the high electronega:vity of Pt, Pt(111) surfaces are nega:vely charged1 .Consequently, the subsurface is posi:vely charged.The charge then should go to zero in the middle of the slab, where bulk condi:ons are achieved.In a 3-layers system, such as the one discussed in Figure 2, the middle layer corresponds to the subsurface to both surfaces.Our 5-layers model is more in line with this picture (Figure 1b, top).It has to be noted that the use of Bader charges introduces some inaccuracies in the assignment of the charge to atoms.Indeed, its magnitude, is the result of the integra:on over atomic volumes which might include tails of the charge actually belonging to neighbour atoms, .PloUng the original average total charge density distribu:on along Z, however, we can see that it goes to zero in the middle of the metal slab (Figure 1, bo@om), even for the three layer system.When ∆ = 0  the total excess Bader charge on both slabs (the sum of the charge on all layers) goes to 0.0 ± 1 $% |e|.When ∆ ≠ 0  the total excess Bader charge becomes nega:ve on one slab and posi:ve one the other slab.The plot of the excess Bader charge per layer reported in Figure 2 served to show that such varia:on happens only on the two internal surfaces, as expected for a capacitor.Since the ini:al charge of the internal surfaces at ∆ = 0  is nega:ve, these layers will only become "more nega:ve" and "more posi:ve" with respect to their ini:al state.It is the total excess Bader charge on the slabs that was used for the calcula:on of the capacitance.
Consequently, the HP-DFT can correctly capture the behaviour of the plate.
To clarify this, we have modified the manuscript as follows: lines 200 to 214 of main_reviewed.pdf:"When ∆ = 0 eV the total excess Bader charge on both slabs (obtained from the sum of the excess Bader charge on each atom) goes to zero (within 0.0 ± 1 $% |e|).When ∆ ≠ 0 eV, the total excess Bader charge becomes negaGve on one slab and posiGve on the other.As can be observed in the figure, such variaGon in charge only occurs on the internal surfaces of the capacitor (3rd and 4th layers) as expected.Due to the high electronegaGvity of Pt, the atoms on the surfaces are negaGvely charged even at ∆ = 0 eV.Therefore, when ∆ ≠ 0 eV the average charge per atom on the internal surfaces will become "more negaGve" or "more posiGve" with respect to the zero bias case.To compensate, the charge of the subsurface (2 nd and 5 th layers in our 3-layer slabs), will have a posiGve charge.Figure (S2e) in the SupporGng Material shows the average Bader charge per atom in each layer of the plates for the Pt(111)(2×2×5) system.It can be observed that in the 5-layer plates, while the subsurfaces (2 nd and 4 th layers for the led plate, and 7th and 9th layers for the right plate) are sGll posiGvely charged, the charge tends zero in the middle of the slab (3 rd and 6 th layers).This observaGon indicates a realisGc descripGon of the metallic slabs.The overall charge on the slabs for the Pt(111)(6×6×3) system is given in Table 2.".
We have also added the following text and

It is worth noGng that the magnitude of the average charge per atom is a result of Bader analysis, which integrates over atomic volumes. Plofng the average total charge density distribuGon along Z, however, we can see that the total charge goes to zero in the middle of the metal slab for all systems, as shown in Figures (S2g), (S2h) and (S2i).".
5 The electrochemical stability window of water is 1.23 V, but the electrochemical poten:al range calculated in the manuscript is -4 eV to 4 eV, far exceeding the electrochemical window of water molecules.At such high volt ages, water molecules should have already undergone reac:ons.
For water molecules to react and split there is a barrier to be overcome.Since this is a geometry relaxa:on at 0K, there is not enough energy to break bonds in molecules.The voltage is enough to cause the molecules to overcome the barrier to rotate to H-down and to even desorb but not to break bonds.

Reviewer #2:
6 The study conducted by Buraschi and colleagues explores the polariza:on effects of the electric double layer at the Pt(111)/water interface.The main novelty of this work is integra:ng a constant poten:al method originally developed by Zauchner and co-workers ((PRB, 2018, 97, 045116)) into CP2K.That method facilitates the examina:on of electrochemical interfaces at predetermined chemical poten:als, a feature of significant interest in fields such as electrocatalysis, corrosion, and electropla:ng.(a) Despite the a@rac:veness of the constant poten:al method for modelling electrochemical interfaces, the manuscript does not clearly delineate the main advantages of the used method over exis:ng constant-poten:al methods, such as those proposed by Bonnet, PRL, 2012, 109, 266101 and Bouzid, JPCL, 2018, 9, 1880.(b) Furthermore, the inves:ga:on into water orienta:on and adsorp:on energy in response to electrode poten:al, u:lizing a simplified water bilayer model, does not convincingly reveal novel electrochemical insights.Consequently, this work appears to be more methodologically focused, sugges:ng its suitability for a journal with a technical emphasis.
(a) See our response to Q1 (b) The water bilayer system did indeed start as a simple model for proof of concept for the methodology.Its study with HP-DFT calcula:ons, however, unveiled relevant insight for electrochemistry.Even more relevant, is that this insight came from such a simple system, without the need for computa:onally expensive model and/ or formalism, highligh:ng the importance of direct poten:al control.
Our study finds that the applied EC poten:al significantly influences coverage, interface structure and adsorp:on, aligning with prior research conducted through complex models and costly computa:onal methods such as AIMD, NEGF, GC-DFT, etc. Therefore HP-DFT emerges as an efficient tool, capable of accurately capturing polariza:on responses with small models, thus refining our understanding of poten:al effects.
7 In this work, the terminology "constant poten:al" is used to describe the poten:al difference between a working electrode and a counter electrode.However, it is not the case of electrochemistry, in which the "constant poten:al" should be the poten:al difference between a working electrode and a reference electrode.I am not convinced how this method can be used for studying electrocatalysis at the so-called constant poten:al condi:on.
We recognise that "constant poten:al" could indeed be a misleading terminology.In the HP-DFT formalism we fix the electrochemical poten:al of the electrodes.We can fix the electrochemical poten:al of two plates to create an electrochemical poten:al difference, which is a constant poten:al difference.This directly models the behaviour of electrochemical capacitors and junc:ons.In electrocatalysis the electrode poten:al is measured experimentally with respect to a reference electrode.To define the electrode poten:al drop in a half-cell -the Galvani poten:al, which although not experimentally measurable, can be evaluated with a calcula:on -it is necessary to define a reference level.Following founda:onal work of TrasaU and Sprik/Chang, this could be represented by the vacuum level in front of the electrolyte solu:on or the bulk electrolyte poten:al, respec:vely.
In our system, the water layer is not sufficiently thick to effec:vely shield the surface dipole or provide a reference bulk electrolyte level, leading to inaccuracies in establishing such a reference level.However, this does not represent a fundamental limita:on of our methodology, preven:ng to use it for studying the electrocatalysis of half cells at constant poten:al (provided that an appropriate reference point is accurately iden:fied within the system).Despite this limita:on, it is remarkable that we s:ll observe behaviours induced by the modifica:on of the local interfacial electric field with bias, that are consistent with what observed with much larger models.More in general, to obtain a correct descrip:on of the poten:al drop, two levels must be iden:fied anyway (not only in our methodology).Our approach allows to flexibly control one of the two, adding precious flexibility to the modelling of interfacial electrochemical behaviours with respect to other methodologies.Interes:ngly, this reference point can be used to place the electrode poten:al on an absolute scale.More specifically, to accurately define a reference level, we plan to use this approach using systems with a full water layer.Further studies are ongoing.However, unravelling the intricacies of the defini:on of an absolute scale for the electrode poten:al is out of the scope of present proof-ofprinciple paper.
Following this comment, we have modified the manuscript, changing the terminology "constant poten:al" to "constant electrochemical poten:al" throughout.See also modifica:ons introduced below.
8 The calcula:on of water adsorp:on energy, as outlined in Equa:on (6), raises ques:ons.Typically, the reference state for water would be its gaseous form, making the inclusion of the last term in Equa:on (6) ques:onable.A revision or jus:fica:on for this term would be beneficial.
The reviewer makes a good point about the equa:on we used.We have modified the formula for the adsorp:on energy as follows: The term explicitly containing the EC poten:al of the slab has been removed.This is Equa0on (8) of main_reviewed.pdf.
The new values of adsorp:on energy resulted in the following plot (Figure (5) of main_reviewed.pdf): The new equa:on represents the energy for inser:on of a water bilayer into the charged capacitor.
The original equa:on stemmed from the idea that  could poten:ally be used as the electrochemical poten:al of the reference electrode, as discussed in Q2(c).Consequently, the manuscript has been modified as follows: lines 249 to 253 and lines 261 to 268 of main_reviewed.pdf:"The adsorpGon energies for the for inserGon of a water bilayer into the charged capacitor, ∆ &'( (∆), were calculated according to: EquaGon (8), where  (,(-./(∆) is the total energy of the capacitor+water system at ∆,  (0&1 (∆) is the energy of the slabs without water at ∆,  )*+ () is the energy of the molecule in gas phase and  )*+ is the number of adsorbed water molecules." We would like to address how this methodology can be used for electrochemical applica:ons.In EC studies the electrode poten:al needs to be explicitly included in the evalua:on of adsorp:on energies.Consequently, a sensible reference level needs to be established.This can be done, for example, in situa:ons where an extended bulk region is included in the simula:on, where the poten:al of the bulk water could serve as the reference level.Such a level, however, does not exist in our systems, as there is not enough water to screen the dipole.Consequently, we have evaluated the adsorp:on energy of the water layer on a half cell, using as a reference the best evalua:on for the vacuum level in our system.For the sake of this discussion, we will take the Hartree poten:al of in the middle of the vacuum region (where the average Hartree poten:al profiles for all the systems cross) as the reference value to calculate the "work func:ons" of the lep and right electrode.Indeed, if we applied a dipolar correc:on in the vacuum region to cancel out the field, the resul:ng vacuum level would correspond to the value of the reference level indicated in figure (4).These "work func:ons" therefore are calculated as:  " =  /2''0.−  " ( 2 ) and  # =  /2''0.−  # ( 3 ) respec:vely, where  /2''0. in the Hartree poten:al of in the middle of the vacuum region and  " and  # are the local Fermi levels on the lep and right plate respec:vely.The inclusion of these terms in the equa:on for the adsorp:on energies then becomes:    !"##$% the surface during geometry op:miza:on.Such reduc:on in water coverage further stabilizes the bilayer making it more stable than its full coverage H-down counterpart.This observa:on aligns with recent literature in showing that capaci:ve response of the interface to a posi:ve EC poten:al is primarily driven by the increase in surface coverage of posi:vely charged water molecules.
In a system with a well-defined and sensible reference level, the same equa:ons (2) to (4) can be applied using such level instead of  /2''0. .The above discussion has been men:oned to the manuscript: lines 284 to 296 of main_reviewed.pdf:"To demonstrate the applicability our methodology to electrochemical studies, we have also evaluated the adsorpGon energy of the water bilayer on a half cell, using as a reference the best evaluaGon for the vacuum level in our system.In this case, the agreement with AIMD studies is even more remarkable, as evidenced by the increased stabilizaGon of a higher number of water molecules at elevated potenGals.This observaGon, in fact, aligns with recent literature in showing that capaciGve response of the interface to a posiGve EC potenGal is primarily driven by the increase in surface coverage of posiGvely charged water molecules . 3,17,50.This discussion is presented in more details in the SupporGng Material." The detailed discussion presented above is also given in the SI, lines 114 to 148 plus FigureS6 of Suppor0ng_reviewed.pdf.
9 To be@er elucidate water polariza:on at charged surfaces, considering a water monomer model might offer a more straighuorward analysis.
The analysis of the effects of the applied EC poten:al on a single molecule adsorbed is showcased in the Suppor:ng Info.
We have modified the SI as follows: • lines 50 to 54 of Suppor0ng_reviewed.pdf:"The more posiGve the right surface becomes the more electronic charge H2O donates to the electronegaGve Pt surface.Conversely, as the led electrode surface becomes more negaGvely charged, less electronic charge is donated to Pt; consequently, the molecule becomes less and less posiGvely charged, going towards zero (the charge of the isolated molecule).18 I have difficul:es with the formula:on "the poten:al increases" (p.11).I assume the authors mean "becomes more posi:ve/nega:ve"?
Yes, reviewer 3 is correct.We have modified the manuscript to make this point clearer.See our reply to Q5 for the modifica:ons to the manuscript.
19 What does a "full coverage water bilayer" mean -a 12 water molecules structure in the first and 12 water molecules in the second layer?What is the meaning of "low water coverage structures"?Please be more precise.
As correctly pointed out by the reviewer, the full coverage water bilayer is the system with 12 molecules on the 1 st and 12 molecules in the 2 nd adsorp:on layer.This system is defined "full coverage" as 12 is the maximum number of molecules which can be chemisorbed on our structures (6x6 surface) for a H-down and H-up bilayer structure.The "low water coverage structure" are those were some of the water molecules were removed from the 1 st adsorp:on layer.
We have modified the manuscript to make this more understandable: lines 131 to 135 of main_reviewed.pdf:"This system is defined as "full coverage" as 12 is the maximum number of molecules which can be chemisorbed on the 6x6 surface for an ice-like H-down and H-up bilayer structure.AddiGonally, we created "low water coverage structures" by removing 2 and 4 water molecules from the 1 st water layer in both H-up and H-down models.These structures are denoted as 12H2O, 10H2O and 8H2O.".
20 I believe that a cita:on is missing on p.12 when the authors say that their "results also align with recent literature".
We 21 What is the meaning of E_h?
E_h refers to Hartree energy, the standard units used in many computa:onal codes; it is the unit of energy in the atomic units system.
22 Please be more specific when wri:ng that the "HP-DFT formalism in calcula:ons did not exhibit a significant decelera:on of the SCF cycle."Some numbers were provided in the SI comparing :mings for one cycle.What is not clear to me is, how long does the overall equilibra:on of the system take.Aper all, one needs to equilibrate the various leads.How efficient the overall performance is, would likely very much depend on the implementa:on and used algorithm.
See our response to Q1 and Q14(c) for addi:onal detail on code efficiency.
In addi:on, we would like to clarify that the leads are not a physical part of the system as in NEGF approaches, so they don't need to be equilibrated (A HP-DFT calcula:ons is done in one step).
23 While I can guess which DOS in Fig. 2b relates to the lep and which to the right lead, it would be helpful to also write this explicitly in the figure.
We have modified the figure according to the reviewer's sugges:on.The distances are measured from the centre of mass of the inner layers of the plates.The distance used for the calcula:on of the capacitance, however, has been corrected for the Van der Walls radius.Following the reviewer's ques:on, we have clarified this point in the.See our response to Reviewer 2, Q11, also repor:ng our modifica:ons to the manuscript.
Figure to the SI: lines 26 to 41 and Figure S2 of Suppor0ng_reviewed.pdf:"Figures (S2a), (S2b) and (S2c) show the average Hartree potenGal along the Z-axis for the Pt(111)(2×2×3), Pt(111)(2×2×5) and Pt(111)(6 ×6×3) systems, respecGvely.The local Fermi levels of the led and right plates are idenGfiable in each model.Figures (S2d) to (S2i) illustrate the charge distribuGon within the plates of each model.Due to the high electronegaGvity of Pt, Pt(111) surfaces are negaGvely charged 1-3 and, as a consequence, the subsurfaces are posiGvely charged.The charge then should go to zero in the middle of the slab, where bulk condiGons are achieved.In our 3-layer systems, Figures (S2d) and (S2f), the middle layers for both slabs (2 nd and 5 th layers) are also the subsurfaces to both surfaces.Consequently, they are posiGvely charged to compensate the negaGve charge on the surface.In our 5-layer model, while the subsurfaces (2 nd and 4 th layers for the led plate, and 7 th and 9 th layers for the right plate) are sGll posiGvely charged, the charge tends zero in the middle of the slab (3 rd and 6 th layers), as shown in Figure (S2e).

Figure 2 .
Figure 2. Adsorp<on energies for the water bilayer plo?ed as a func<on of the poten<al at different coverages: we compare both H-up and H-down configura<ons for cases with 12H2O, 10 H2O, and 8 H2O chemisorbed water molecules.
)*+ { (,(-./(∆) −  (0&1 (∆) −  )*+ *  )*+ () + [ " (∆) *  " (∆) +  # (∆) *  # (∆)] } ( 4 )where  " (∆) and  # (∆) are the "work func:ons" of the lep and right electrode at ∆ and  " (∆) and  # (∆) are the excess Bader charges on the lep and right plate at ∆.The resul:ng adsorp:on energies are reported in the following plot:The effects of including the term between square brackets in Equa:on (4) are presented in Figure(3).The data in Figure(3) indicates that the low coverage H-down structures are s:ll the most stable ones for values of ∆ between −1  and +1 .At ∆ = 4 , however, the full coverage Hdown water bilayer becomes the most stable structure.At ∆ = −4 , instead, the full coverage H-up bilayer structure is the most stable one.In this system, most of the molecules in the 2 nd adsorp:on layer flipped from an H-up to an H-down configura:on, stabilizing the bilayer.This is also the structure where one of the chemisorbed molecules in the 1 st layer spontaneously desorbs from

Figure 4 .
Figure 4. Adsorp<on energies for the water bilayer calculated with the equa<on (4).
". • lines 71 to 77 of Suppor0ng_reviewed.pdf:"Figure (S4) shows the electron density difference between the system at ∆ = 4  and the system at ∆ = 0 .An accumulaGon of electronic charge (indicated by the blue isosurface) between the water and the led plate, suggesGng a stronger bond between the molecule and the surface.Conversely, a depleGon of electronic charge (indicated by the orange isosurface) between the water and the right plate indicates a weakening of the molecule-surface bond.Both observaGons are in agreement with the data obtained from the Bader and geometry analysis of the system."• We have also added the following picture, which shows how the charge polarises under applied EC poten:al difference, as Figure S4 in the Suppor:ng Info: As can be observed from Figure (6), the overall trend remains the same as that of the ∆ &'( (Figure(2)).Therefore, our discussion on the structure and stability of the bilayer remains unaffected.

Figure 6 .
Figure 6.Adsorp<on energies for the water bilayer plo?ed as a func<on of the poten<al at different coverages: we compare both H-up and H-down configura<ons for cases with 12H2O, 10 H2O, and 8 H2O chemisorbed water molecules.
Figure (3b) of main_reviewed.pdf.24 It would be helpful to incorporate the labels of the graphs in Fig. 4b into the figure.We have modified the Figure according to the reviewer's sugges:on.Figure (6b) of main_reviewed.pdf.25 How are distances measured?Do the authors use the centre of mass of a water molecule or maybe the oxygen in a water molecule?

Figure S1, which
compares number of SCF cycles and average :me per SCF cycle steps in standard DFT and HP-DFT calcula:ons, has been added to Suppor0ng_reviewed.pdf.